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Abstract 

The bifurcation diagram for a vibro-fluidized granular gas in connected compartments is 
constructed and discussed. At vigorous driving, the uniform distribution (in which the gas is equi- 
partitioned over the compartments) is stable. But when the driving intensity is decreased this 
uniform distribution becomes unstable and gives way to a clustered state. For the simplest case, 
N = 2, this transition takes place via a pitchfork bifurcation but for all > 2 the transition in- 
volves saddle-node bifurcations. The associated hysteresis becomes more and more pronounced for 
growing A^. In the bifurcation diagram, apart from the uniform and the one-peaked distributions, 
also a number of multi-peaked solutions occur. These are transient states. Their physical relevance 
is discussed in the context of a stability analysis. 

PACS numbers: 45.70.-n, 02.30.Oz 
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I. INTRODUCTION 



One of the key features of a granular gas is the tendency to spontaneously separate into 
dense and dilute regions [|I], 0, ||, H, |^. This clustering phenomenon manifests itself in a 
particularly clear manner in a box that is divided in a series of connected compartments, 
with a hole (at a certain height) in the wall between each two adjacent compartments. The 
system is vibro-fluidized by shaking the box vertically. With vigorous shaking the granular 
material is observed to be distributed uniformly over the compartments as in any ordinary 
molecular gas. Below a certain driving level however, the particles cluster in a small subset 
of the compartments, emptying all the others. 

For N = 2 the transition from the uniform to the clustered state is of second order, taking 
place through a pitchfork bifurcation . For = 3 it was recently found that the transition 
is hysteretic. It is a first order phase transition, involving saddle-node bifurcations P] . This 
difference has been explained by a flux model. In the present paper we will use the same 
flux model to construct the bifurcation diagrams for arbitrary A^. 

The main ingredient of this model is a flux function F{n), which gives the outflow from a 
compartment to one of its neighbors as a function of the fraction of particles (n) contained 
in the compartment 0. The function F{n) starts out from zero at n = and initially 
increases with n. At large values of n it decreases again because the particles lose energy 
in the non-elastic collisions, which become more and more frequent with increasing particle 
density. So F{n) is non-monotonic, and that is why the flux from a well- filled compartment 
can balance that from a nearly empty compartment. 

Assuming that the granular gas in each compartment is in thermal equilibrium at any 
time (in the sense of the granular temperature [§]) the following approximate form for F{n) 
can be derived [^: 



which is a one-humped function, possessing the features discussed before (See Fig. ||). In 
the above equation is the fraction of particles in the k-th compartment, normalized to 
^ nfc = 1. The factors A and B depend on the number of particles and their properties (such 
as the radius, and the restitution coefficient of the interparticle collisions), on the geometry 
of the system (such as the placement and form of the aperture between the compartments). 




(1) 
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and on the driving parameters (frequency and amplitude). The factor A determines the 
absolute rate of the flux, and will be incorporated in the time scale, which thus becomes 
dimensionless. The clustering transition is governed only by B. 

The time rate of change fik of the particle fraction in the A;-th compartment is given by 
the inflow from its two neighbours minus the outflow from the compartment itself, 

hk = F{nk-i) - 2F{nk) + F^+i), (2) 

with k = 1,2, ..,N. Here we have assumed that the interaction is restricted to neighboring 
compartments only. 

For a cyclic arrangement the above equation is valid for all N compartments (with k = 
+ 1 equal to = 1). If we take non-cyclic boundary conditions, by obstructing the flux 

between two of the compartments, the equation has to be modified accordingly for these 

compartments. 

The total number of particles in the system is conserved (^^n^ = 1), so 

J2^k = 0. (3) 

k 

Statistical fluctuations in the system would add a noise term to Eq. (^, but we will not 
consider such a term here. So the present analysis has to be interpreted as a mean field 
theory for the system. 



Equation 



can also be written in matrix-form, as n = M ■ F, or more explicitly: 
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The given matrix M corresponds to a cyclic arrangement of the compartments. A 
similar matrix can be written down for the case of a non-cyclic arrangement. We will come 
back to this later, when we will see that most of the results for the cyclic arrangement carry 
over to the non-cyclic case. 

It is easily seen, from the fact that the elements of each row of M sum up to zero, that 



1) is an eigenvector. The corresponding eigenvalue A = physically reflects 
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FIG. 1: The solutions n_ and n+ of F{nk) = constant, cf. Eq. |5[ Also shown is how the flux 
balance responds to an increase of n_ by an amount of 6n (see also Eq. |Tl|). The diagram on the 
right hand side depicts the relation between F and the quantity a = F'{n^)/F'{n^), which plays 
an important role in the stability analysis of section 3. 



the fact that the compartments cannot all be filled (or emptied) simultaneously: ^j^, = 
or Mik = 0. For future reference we note that all the other eigenvalues of M are negative 
(see Appendix). 

The remainder of the paper is set up as follows. In Section II we show how to construct 
the bifurcation diagram, on the basis of Eq. (^, for an arbitrary number of compartments. 
In Section III we discuss the stability of the various branches in the diagram. Section IV 
discusses the physical consequences resulting from the diagram, in particular in the limit 
for N oo. Finally, Section V contains concluding remarks. The paper is accompanied by 
a mathematical Appendix, in which some essential results concerning the stability analysis 
are derived. 



II. CONSTRUCTING THE BIFURCATION DIAGRAM 

To calculate the bifurcation diagram, we have to find the fixed points of Eq. (^) as a 
function of the parameter B, i.e., those points for which ri^ = M ■ F = 0. So F must be a 
multiple of the zero-eigenvalue vector 1 = (1, 1, ... , 1). This tells us that, in a stationary 
situation, all components of the flux vector F are equal: there is a detailed balance between 
all pairs of neighboring compartments. This rules out, for instance, the possibility of stable 
standing-wave-like patterns with equal but non-zero net fluxes throughout the system. The 
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fixed point condition now becomes 

FinA = constant 

(5) 

Since F is a one-humped function, F{nk) = constant fias two solutions, wliicli will be 
called n_ and n+ (see Fig. |l|). Every fixed point can be represented as a vector with elements 
n__ and n^ (in any order, and summing up to 1) corresponding to a row of nearly empty 
and well-filled compartments. Let us call the number of well-filled compartments m. Apart 
from the ordering of the elements, every fixed point is then specified by only two numbers: 
n+ and m. 

Before actually calculating the bifurcation diagram, it is convenient to replace the fraction 
n by the (also dimensionless) variable z = Nny/B, as then the flux (^ simplifies to F{zk) oc 
z^exp(— 2;^). The fixed point condition Eq. (|^) then reads: 

Fizk) = constant 

(6) 

So the -B-dependence has been transferred from F to the particle conservation, and this 
enables us to determine the entire bifurcation diagram from one single graph. This is 
illustrated in Fig. Q for the case of A^=5 compartments. 

First, the one-humped function F{z) is inverted separately on both sides of the maximum, 
yielding the functions Z-{F) and z^{F). Then, we construct the sumfunctions: 

Sm{F)=mz4F) + {N-m)z4F) (7) 



Now, from Eq. (^), the fixed points are found by intersecting the horizontal line z = NyB 
with the sumfunctions Sm{F). In Fig. ^ this is done for S=1.08. Each intersection point 
yields a pair or equivalently {n_,n^}. Repeating the procedure for all B, we 

obtain the bifurcation diagram depicted in Fig. |^. 

It contains several branches. First, a horizontal line (from the sumfunctions Sq and S'5) 
corresponding to the equal distribution n+ = n_ = 0.2 = Second, the branches 

corresponding to the m = 1 clustered state (from Si), which aX B = 1 goes over into the 
m = 4 state (from 5*4). And third, the branch of the m = 2 clustered state (from 6*2), which 
at i? = 1 becomes the m = 3 state. The physical appearance of these solutions is sketched 



5 




FIG. 2: Inverted flux functions z^{F) and z^{F) and the + 1 sumfunctions SmiF) = mz^{F) + 
(A^ — m)z-.[F),m = 0,1,... ,N. Here we picked N = 5. The points of intersection with the 
horizontal line z = N\/^ represent the fixed points for the parameter value B. Curves 5*0 and S5 
correspond to the uniform distribution (below and above the critical point B = 1, respectively) 
and the other curves belong to clustered states. Note that Sq joins smoothly with S5 at i? = 1 (i.e. 
z = Ny/B = 5), and so does 5i with 5*4, and 5*2 with 5*3. 



in the small diagrams. Note that only the m = branch (i.e., the uniform solution up to 
B = 1) and the outer m = 1 branch are stable. All the other branches are unstable, as will 
be discussed in the next section. 

At B = 1, where the branches intersect with the uniform distribution = n„ = 
we have a critical point. In the flux function one passes the maximum here. This means 
that n+ and n_ are switched (relatively empty compartments become relatively filled, and 
vice versa), so m-branches change into (A^ — m)-branches. From a physical point of view, 
the most important thing that happens at the B = 1 intersection point is the destabilization 
of the uniform distribution. 

The saddle-node bifurcations of the m = 1 and m = 2 branches correspond to the minima 
of the sumfunctions 5*1 and S2 respectively, which in Fig|^ can be seen to occur at F ^ 0.014 



for 5*1 and F ^ 0.202 for 5*2. In general, if a sumfunction Sm{F) has a minimum for a certain 
B, the associated m branch will have a bifurcation. So the bifurcation condition is that the 
derivative dSm{F) / dF equals zero, or equivalently: 



dz+ \ N — m 



dF 



(8) 



Not surprisingly, the quantity on the left hand side {dz_./dz^ = a) will play an important 
role in the stability analysis of the next section. 



III. STABILITY OF THE BRANCHES 



The stability of the branches (i.e., of the fixed points) is determined by the eigenvalues 
of the Jacobi matrix J corresponding to Eq. (^), with components: 

^.-^=Ew.)|:^=M..f'K) (9) 

Here F' denotes the derivative of F with respect to n. Note that the Jacobi matrix can also 
be written as the product of M and the diagonal matrix D = diag(-F'(r;,i), . . . ,F'{n]sf)), 
see also Eq. ( |T8|) in the Appendix. For a fixed point the only diagonal elements that occur 
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FIG. 3: Bifurcation diagram for = 5. It has been obtained from Fig. || by converting, for all B, 
each z+}-pair belonging to a point of intersection to a {n_, n+l-pair. Note that all branches 
come together at the critical point B = 1. The (stable) m = branch becomes the (unstable) 
m = 5 branch, the m = 1 branch turns into the m = 4 branch, and m = 2 switches to m = 3. 
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are F'{n+) (m times) and F'{nJ) {N — m times), in any order. The ratio between these 
two functions is precisely the quantity we encountered earher in the bifurcation condition 
Eq. (H), namely cr: 

The Jacobi matrix J has eigenvalues, one of which is always zero. The other — 1 
eigenvalues depend on m and the value of a. 

For m = (the equipartitioned state) all non-trivial eigenvalues are negative, up to the 
point B = 1. This can be seen either by direct numerical calculation, or analytically (see 
Appendix). At i? = 1, the m = state becomes the m = N state. Here, the functions 
F' in the Jacobi matrix (^ change sign, and so do all of its eigenvalues. So suddenly the 
uniform state has A^ — 1 positive eigenvalues, which implies a high degree of instability. 
Only in the limit B ^ oo does the uniform state regain some of the lost terrain: the 
magnitude of all positive eigenvalues tends to zero here. Physically speaking, in this limit 
the vibro-fluidization is too weak to drive the particles out of the boxes anymore. 

As for the other values of m, in Fig. |^ we have plotted the numerically evaluated eigen- 
values (as functions of a) for the system with A^ = 5 compartments. 

For m = 1, we see that there are three eigenvalues that are always negative. The fourth 
non-trivial eigenvalue changes sign at cr = —0.25. This corresponds to the saddle-node 
bifurcation of the m = 1 branch in the bifurcation diagram (Fig. ^), and the bifurcation 
value of a is in agreement with Eq. (^. The region to the right of a = —0.25 (where 
all non-trivial eigenvalues are negative) belongs to the stable outer branch. The left part 
a < —0.25 belongs to the unstable inner branch, up to the point 5 = 1 (at o" = — 1), where 

the m = 1 branch goes over into the m = 4 branch. That is, the state {H } now 

switches to { — 1- + + +}. At the same time all eigenvalues change sign, so suddenly we have 
3 positive eigenvalues, which is only one less than for the uniform m = 5 state. (Indeed, the 
only stable manifold of the m = 4 fixed point comes from the direction of the completely 
unstable m = 5 state). The positive eigenvalues never cross zero anymore (there are no 
bifurcations beyond B = 1) but, as before, in the limit S — oo (o" — 0) they go to zero. 

For m = 2 there are two possible configurations: {+ H } and {H 1 }. Due 

to the cyclic symmetry, all other combinations are equivalent to these two. The eigenvalues 
of the first configuration are given by the dotted lines, and those of the second by the solid 
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m = I 



m = 2 




FIG. 4: Eigenvalues of the Jacobi matrix J as a function of a, for the branches m = 1, 2, 3, 
and 4. Rather than plotting Aj, we display Aj = Xi/ F' (n^), because this yields a more clear-cut 
picture. Negative eigenvalues represent stable directions of the branches, and positive eigenvalues 
represent unstable ones. A zero crossing (such as for m = 1 and m = 2) indicates the occurrence 
of a saddle- node bifurcation. The value cr = corresponds to the limit B ^ oo, and o" = — 1 to the 
critical point B = 1. At this point, the eigenvalues of m = 1 and m = 4 are equal but opposite in 
sign: the transition from the one branch to the other is marked by a distinct drop in stability. The 
same is true for the eigenvalues of m = 2 and m = 3, and also (not depicted) for those of m = and 
m = 5. For m = 2, 3 there are two different cluster-configurations, with different eigenvalues. The 

dashed lines correspond to {+ H } for m = 2, which goes over into { h + +} for m = 3. 

The bold lines apply to the slightly more stable configurations {H 1 } and { — I h +}. 



lines. Although they are very similar (and are represented by exactly the same branch in 
the bifurcation diagram), it is clear that the second configuration is the more stable of the 
two. Apparently the two well-filled compartments prefer to keep a distance. 

The saddle- node bifurcation of the m = 2 branch takes place at cr = —2/3 [cf. Eq. (|^)], 
where the third non-trivial eigenvalue goes through zero. The fourth non-trivial eigenvalue 
always remains positive, indicating that the m = 2 branch never becomes completely stable. 
(As a matter of fact, only the m = branch and part of the m = 1 branch can be completely 
stable). Note that for o" ^ (large B) the positive eigenvalue tends to zero, so the degree 
of instability is quite weak there. 
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FIG. 5: Bifurcation diagram for = 6. Note the pitchfork bifurcation at i? = 1. 



At B = 1 the m = 2 branch becomes the m = 3 branch, with the two configurations 

{ h + +} and { — I h +}, and with all eigenvalues switching sign. As we see, the more 

dispersed configuration is again the less unstable one. Also the phenomenon of all positive 
eigenvalues going to zero as a approaches zero (the weak driving limit i? ^ oo) is again 
apparent. 

In the present example for N = 5, and in fact for all odd values of A^, the branches in the 
bifurcation diagram are all born by means of a saddle-node bifurcation. But for even values 
of A^ this is different: in that case there is one branch-pair that springs from the uniform 
distribution, at i? = 1, by a pitchfork bifurcation. This is illustrated in Fig. |^ for A^ = 6. 
Here one sees all the branches that were present already for A^ = 5, only slightly shifted 
towards the left, plus an additional pair of branches (m = 3) bifurcating in the forward 
direction from B = 1. 

The special status of the branch m = N/2 is also evident from Eq. (H), which tells us 
that the bifurcation condition for this branch is a = — 1. This condition is fulfilled only by 
n+ = ?2_ = 1/^/B = 1/N. So, unlike all other branches, this one originates at B = 1 from 
the (until then stable) uniform state. Related to this, the branch is the only one that is 
symmetric for interchanging n+ and n_. 



IV. PHYSICAL ASPECTS 



The bifurcation analysis from the previous section can also be understood from a more 
physical point of view. To this end, let us first have a closer look at a 2-box system. In the 



10 



equilibrium situation the net flux between the two boxes is zero, with one fllled (n+) and 
one nearly empty box. Suppose the level of the empty box is raised by an amount 5n. 

The level of the fllled box then decreases by an equal amount and the net flux ,_|_ from 

the empty to the fllled box becomes (see also Fig. |I]): 

= F(n_ + 5n) — — 5n) 

dF dF\ ^ ,^ ^ dF ^ (11) 
+ -; — ] Sn = (1 + a) - — 6n 



du- dn+ J dn. 

where we have used that a = dn^/dn^ and neglected the higher order terms in the Taylor 
expansion. There are two different regimes. If a > — 1, the net flux is positive (as F'inJ) is 
always positive), so particles are flowing from n_ to ?7,+ , restoring the equilibrium position. 
This is actually the situation along the entire m = 1 branch, for all 1 < 5 < oo. For cr < — 1 
(a situation which does not occur for our choice of F), the net flux would be negative, 
raising the level of the emptier box even further, away from the equilibrium position. In the 
borderline case, a = —1 (at B = 1), the system is indifferent to inflnitesimal changes. 

This argument is readily generalized to the A^-compartment system, for an equilibrium 
with m fllled boxes. Now we raise the level of all — m nearly empty boxes simultaneously 
by 6n. This is done by lowering all levels in the m fllled boxes by an equal amount, which 
by particle conservation must be equal to 6n{N — m)/m. The equivalent of Eq. ([TT|) for the 
flux between any of the empty boxes to a neighbouring fllled box then reads: 

A'^ — m 

= F(n^ + 6n) — F(n^ 6n) 

m 

dF N -m dF 



- dn (12) 



dn_ m dn^ 

m \ dF 

cr - — On 



N — m J dri- 

From this expression it follows that the transition between a (relatively) stable [a > 
~m/{N — m)) and a (relatively) unstable {a < —m/{N — m)) conflguration is marked 
by the bifurcation condition Eq. (^. So, by straightforward physical reasoning we have 
reproduced the exact result obtained earlier from an eigenvalue analysis. 

The pitchfork bifurcation discussed at the end of Section III is especially important for 
N = 2. In this case it is the only non-uniform branch. To be speciflc, it is a stable 
m = 1 branch. This N = 2 case is the only one without any saddle-node bifurcations, 
and consequently it is the only case where the change from the uniform to the clustered 
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situation takes place via a second order phase transition without any hysteresis. For all 
N > 2 the transition is of first order and shows a hysteretic effect that becomes more 
pronounced for growing N. 

In the limit N oo the hysteresis is maximal: the first saddle-node bifurcation takes 
place immediately after B = 0, and this means that there exists a stable m = 1 solution 
over the entire range B > 0. So, if one starts out from this solution (at a certain value of 
B) and then gradually turns down B, one will never witness the transition to the uniform 
distribution. Vice versa, also the transition from the uniform solution to the m = 1 state 
will not occur in practice, even though the uniform distribution becomes unstable at -B = 1. 
If one starts out from the uniform solution (at a certain value of B below 1) and increases 
B, one will witness the transition to a clustered state, but in practice this will always be 
one with a number of peaks. That is, the system gets stuck in a transient state with m > 1, 
even though such a state is not stable (it has one or more positive eigenvalues). 

The fact is that its lifetime may be exceedingly large, since the flux in the neighborhood 
of a peak and its adjacent boxes (which are practically empty) is very small. Furthermore, 
the communication between the peaks is so poor that usually (even for moderate values of 
A^) the dynamics comes to a standstill in a state with peaks of unequal height. 

Another point we would like to address is that practically the transition to a clustered 
state will take place already before B = 1, because the solution is kicked out of its basin 
of attraction by the statistical fluctuations in the system An example is shown in 
Fig. p. Here we see a snapshot for the cyclic system with A^ = 80 compartments, which 
were originally filled almost uniformly, at B = 0.90. The small random fluctuations in the 
initial condition are sufficient to break away from the (still stable) uniform distribution, and 
one witnesses the formation of a number of isolated clusters. In the further evolution these 
clusters deplete the neighbouring compartments and indeed the whole intermediate regions. 
But the peaks themselves, once they are well-developed, do not easily break down anymore. 

V. CONCLUDING REMARKS 

In this paper we have constructed the bifurcation diagram for a vibro-fluidized granular 
gas in A^ connected compartments. Let us now comment upon the result. 

Starting out from B = 0, i.e, vigorous shaking, the equi-partitioned state is for some time 
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FIG. 6: Results from a numer ical solution of Eq. (|) for N = 80, at B = 0.90. Snapshots are taken 
after lO'^, lO'^, 10'* and 10^ timesteps (iterations). Between 100 and 10,000 iterations a clustering 
pattern is seen take shape. Although strictly speaking this is a transient state, the system gets 
stuck in it. 

the only (and stable) fixed point of the system. For increasing B we first come upon the 
m = 1 bifurcation, where the single-cluster state is born. For all > 2 this happens by 
means of a saddle-node bifurcation, creating one completely stable state and one unstable 
state (with 1 positive eigenvalue). The one with the largest difference between n+ and n_ 
is the stabler one of the two states. Strictly speaking, there are N equivalent single-cluster 
states, since the cluster can be in any of the compartments. 

For further growing B we come across the m = 2 bifurcation, where two unstable 2- 
peaked states are created. The state with the largest difference between ra+ and has 1 
positive eigenvalue, and the other one 2. The two peaks can be distributed in (^) ways 
over the N compartments, but as we have seen they are not all equivalent. When the peaks 
are situated next to each other we have a more unstable situation (the positive eigenvalues 
are larger in magnitude) than when the peaks are further apart. This is generally true for 
m-peaked solutions: of the ways in which m peaks can be distributed, the ones in which 
the peaks are next to each other are the least favorable of all. 

For increasing B we encounter more and more bifurcations, where unstable m-clustered 
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states come into existence (each with 1 more positive eigenvalue than the previous one), and 
for large the bifurcation diagram is covered by a dense web of branches. In Fig. |^ this 
is shown for N=80. The last saddle-node bifurcation takes place shortly before B = 1 and, 
for this even value of A^, is followed by a final pitchfork bifurcation (creating the m = N/2 
branch) at B = 1. 




FIG. 7: Bifurcation diagram for N = 80. The hysteresis extends almost all the way down to 
-6 = 0, and there are numerous transient states (cf. Fig. ^). The only strictly stable branches are 
the m = branch (up to B = 1) at Uk = 1/-/V, and the outer m = 1 branches. Naturally, the 
upper m = 1 branch approaches n^. = 1, the upper m = 2 branch approaches n/^ = 1/2, the upper 
m = 3 branch = 1/3, etc. The overlay picture shows the neighborhood of the critical point at 
B = l,nk = 1/N in more detail. 

The uniform solution (or m = state) is stable until B = 1, with — 1 negative 
eigenvalues and 1 zero. For i? > 1 all its negative eigenvalues become positive, making it 
suddenly the most unstable state of all. Also, it now formally becomes the m = N state. 
Moving away from this uniform solution one encounters first the m = A^ — 1 branch with 
N — 2 positive eigenvalues, then the m = N — 2 branch with A^ — 3 positive eigenvalues, etc. 
Finally, one arrives at the outermost m = 1 branch, which has no positive eigenvalues. This 
is the only solution that is completely stable for i? > 1. But as we have seen in the previous 
section, on its way from the uniform distribution to the single peaked state, the system can 
easily get stuck in one of the transient states (especially for large A^) even though these are 
not strictly stable. 
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Throughout the paper, we have concentrated on the case where the compartments are 
arranged in a cychc manner. But in doing so, we have in fact also solved the non-cyclic case. 
Here we close the hole in the wall between the 1st and A^th compartment, and consequently 
the flux between them is zero. The matrix M then takes the following form [differing from 
the cyclic one only in the first and last row, cf. Eq. (H)]: 



/-II ■■■ \ 
/ 1 -2 1 ■■■ n n \ 



]V[("c) 



1-2 1 








v 








-2 1 

1 -1 



(13) 



The eigenvalue problem for this matrix is treated in the Appendix. One eigenvalue is 
identically zero, and the other A^ — 1 eigenvalues are negative, just like for the cyclic system. 
This leads to a bifurcation diagram that is indistinguishable from that of the cyclic case. 
Even the stability along the branches is the same; only the magnitude (not the sign) of the 
eigenvalues of the Jacobi-matrix J is slightly different for the two cases. 



Finally, it should be emphasized that the results of the present paper do not depend on 
the precise form of the flux function. We have concentrated on the form given by Eq. (P, 
but virtually everything remains true for other choices of this function, as long as it is a 
non-negative, one-humped function, starting out from zero at n = (no flux if there are no 
particles) and going down to zero again for very many particles (no flux also in this limit, 
since - due to the inelastic collisions - the particles form an inactive cluster, unable to reach 
the hole in the wall anymore). Any function with these properties will produce a bifurcation 
diagram similar to that of Eq. (|l]) . 

In the likely case that the range of a = dn^/dn^ is the same, extending from —1 (this 
value is attained in the maximum) to zero (in the outer regions of the flux function, for 
n^\fB 0, n^\fB oo), the bifurcation diagram will have the same number of saddle- 
node bifurcations and the same number of branches. The only things that change are the 
exact position of the bifurcation points, and the magnitude of the eigenvalues along the 
branches. 

Slight differences in the diagram would occur if the slope of F on the n+ side was to 
become steeper than on the n_ side. In that case, the bifurcation condition Eq. (||) would 
also have solutions for m > N/2, thus allowing saddle-node bifurcations for branches with 
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m > N/2. These branches, however, would certainly be quite unstable. 



APPENDIX: ON THE EIGENVALUES OF M AND J 

In this appendix we present the analytical eigenvalues of the flux matrix M [introduced 
in Eqs. and (|13|)] and discuss the eigenvalue problem for the Jacobian matrix J [see 
Eq. (^], thereby determining the stability of the branches in the bifurcation diagram. 

First, we briefly treat the eigenvalues of M. After that, we turn to J. In Subsection 2 
we discuss its zero-eigenvalues: one eigenvalue is identically zero and, by pinpointing the 
zero-crossing of a second eigenvalue, we reproduce the bifurcation condition Eq. (^). In 
Subsection 3 we determine the number of negative eigenvalues of J in the low-driving limit 
(7^0. Likewise, in Subsection 4 we determine the number of positive eigenvalues in the 
(mathematical) limit a — oo. Combining these two results, in Subsection 5, we flnally flnd 
the number of positive eigenvalues of J for general values of a, and this gives the stability 
of the branches over the entire bifurcation diagram. 



1. Eigenvalues of matrix M 

The matrix M in Eq. (^) is closely related to the NxN tridiagonal matrix tridiag(l, — 2, 1) 
associated with the second difference operator known from numerical schemes for solving 
second order pde's. Its eigenvalue problem can be solved exactly |]10[, and the same is true 
for M. The eigenvalues of M are given by: 

A. = -4sin^f^l (14) 



where k runs from to N/2 for even, and from to (A^ — l)/2 for A^ odd. The corre- 
sponding eigenvectors are: 

„.(,)^C.co.((Bi±il^^)+C.=i„((^) (15) 

with i = 1, . . . , N and arbitrary coefficients Ci and C2. 

As we see, the first eigenvalue {k = 0) is zero and the corresponding eigenvector is 
1 = (1, 1, . . . ,1). Physically, this eigenvector represents simultaneous filling of all A^ com- 
partments, and the eigenvalue expresses the fact that this is prohibited (because the 
number of particles in our system is conserved). 
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All non-zero eigenvalues are negative and (except the one for k = N/2 in the case of 
even A^) doubly degenerate. This means that the corresponding eigenvectors span a two- 
dimensional subspace, reflected by the two terms Ci and C2 in Eq. (|1^). Since M is sym- 
metric, and therefore normal, linear subspaces corresponding to different eigenvalues are 
orthogonal. Especially, the eigenvectors of all non-zero eigenvalues span a — 1 dimen- 
sional subspace perpendicular to 1 = (1,1,... ,1). 

The matrix M*^"'^^ for the non-cyclic case, given by Eq. (piSf ), has a different set of eigen- 
values: 

Ar^ = -4sin^(^) (16) 

Here k runs from to A^ — 1. The corresponding eigenvectors are: 

arW = cos(Pi±ilil) (17) 

Just like in the cyclic case, the first eigenvalue equals zero, and all the others are negative. 
However, they are non-degenerate and the corresponding eigenspaces are one-dimensional. 



2. Zero-eigenvalues of matrix J 

Now we turn to the Jacobian matrices. We consider the cyclic version J, with components 
as given in Eq. (H), but the results are also valid for the non-cyclic version. This matrix can 



be written as the product of M and a diagonal matrix D = diag(F'(ni), F'{n2) 



M D 



/ -2F'{nx) F'{n2) 
' F'(ni) -2F'(n2) F'{ns) 
F'(n2) -2F'(n3) 



F'iriN) \ 
» 



(18) 



-2F'{nM-i) F'iriN) , 
F'iriN-i) ~2F'{nM)/ 



\ F'(ni) 

In the context of the bifurcation diagram, the main thing one wants to know is the number 
of positive eigenvalues of J for each branch. This is what we are going to determine now. 

First we note that the eigenvalues of J are real, even though the matrix is not symmetric. 
This is a consequence of the following similarity relationship between J and J^: 



= (M D)^ = D M = D (M D) = D J 



(19) 



This implies that J and have the same eigenvalues, and hence they must be real. 
Because M is singular, J must be too (it has a zero eigenvalue) and so its determinant 
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det(J) is zero. More explicitly: 

det(J) = det(M) ■ det(D) = ^11^^'^)^ det(M) = (20) 
where, for a fixed point with m filled compartments, the product term equals 

For the other eigenvalues we have to look at the characteristic equation det(J — AI) = 0. 
This is a polynomial expression in A, of which the constant term is zero since it is equal to 
det(J). The coefficient L of the linear term is: 

L = J2det{j('''^) = J2(U^) (21) 

k k \l^k J 

where the matrix jC^''^') is the (A^ — 1) x (A^ — 1) matrix obtained from J by deleting its 
/c-th row and its fc-th column. In the right-hand side of this equation, the only product that 
survives is the one that does not contain the trivial (zero) eigenvalue. So: 

L= n (22) 

all non-trivialAj 

Alternatively, the determinant of J*^'^''^) in Eq. ( pi]) can be written in terms of det(M*^^''^^), 
by deleting the /c-th factor from the product in Eq. (pO|): 

L = Y, ( n^'(^') ) det(M('=''=)) (23) 

k \lj^k / 

It can be shown that for all k the determinant det(M('^'^)) is a constant, C, which equals 
{N — 1)(— 1)^^^ in the cyclic, and (—1)^^^ in the non-cyclic case. Thus, Eq. ( P^D reduces 
to: 

k \lj^k / 

For a fixed point with m filled compartments, we can write (using that in the above 
summation each of the products misses either an F'(n+) or an F'{nJ)): 

L = C[F'(r2+)]('"-^)[F'(r2_)](^-'"-^) x 

({N - m)F'{n+) - mF'(n_) j (25) 
= C[F'(n+)]('"-i)[F'(n_)](^-'") ({N - m)a - 



18 



From this equation we conclude that L becomes zero at o" = — m/ (A^ — m). This is exactly 
the bifurcation condition already given in the main text [Eq. (^]. Also, with Eq. (p2D, we 
see that an eigenvalue crosses zero at this value of a. 

It can be shown, by a similar analysis, that the coefficient of the quadratic term is not 
equal to zero at a = —m/{N — m), so not more than one of the eigenvalues changes sign at 
the bifurcation. 



3. Number of negative eigenvalues of J for a ^ 

We now come to the next step in determining the number of positive eigenvalues. We 
again use the definition of a to write: J = F'(r2_)M ■ D, where D = diag(l, ... , 1, a, . . . ,a). 
The factors 1 correspond to the N — m nearly empty boxes and the factors a to the m filled 
boxes. The precise ordering of the factors is not essential for the following argument, so we 
may choose the above order for notational convenience. 

The factor F'{nJ) is always positive, so we only have to deal with M ■ D. Note that only 



D depends on a and that in the limit a ^ this matrix becomes[ll 



lim D = diag(l, . . . , 1, 0, . . . , 0) = P (26) 

P is a projection matrix which projects to the subspace spanned by the first N — m 
unit vectors. It is obviously non-singular, symmetric, and applying it twice gives the same 
result as once: P^ = P. 

Instead of taking the matrix Jq = M ■ P as input for solving our eigenvalue problem (in 
the limit a ^ 0), we will rather look at the matrix P ■ M ■ P which is symmetric and has 
the same eigenvalues as Jq. 

For proof of the last statement, let /i be a (non-zero) eigenvalue of Jq: Jq - x = /ix. Then: 
(P ■ M • P) • (P ■ x) = P ■ (M • P • x) = /i(P ■ x). Note that P ■ x 7^ 0, because otherwise 
also Jq ■ X = M ■ P ■ X would be zero, contradicting the assumption that /i is non-zero. This 
completes the proof. 

The matrix M is negative semi-definite. This means that M has only negative or zero 
eigenvalues or, equivalently, the inner product (x, M ■ x) < for all x. This means that also 
P ■ M • P is negative semi-definite, because: 

(x, P ■ M ■ P ■ x) = (P ■ X, M ■ (P ■ x)) = (y, M ■ y) < (27) 
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In conclusion, Jq has negative and zero eigenvalues only. 

The remaining task is to identify the number of negative eigenvalues, or otherwise stated, 
the rank of the matrix Jq. The statement which we shall prove is that rank(Jo) = rank(P) = 
N — m. 

Proof: Note that the image Im(P) of P is spanned by the first m unit vectors of M^. Its 
kernel Ker(P) is spanned by the remaining N — m unit vectors. Since the kernel of M is 
spanned by the vector 1, the following identities hold: 

Ker(P) n Ker(M) = (28a) 

Im(P) n Ker(M) = (28b) 

Now, for all x G Ker(P) it holds that Jq • x = M ■ (P ■ x) = 0, so Ker(P) C Ker(Jo). 
On the other hand, for all y ^ Ker(P) one has P ■ y = z 7^ 0, with z G Im(P), and 
therefore Jq • y = M ■ z 7^ because of Eq. (|28 b| ) . This means that y ^ Ker(Jo), and thus 
Ker(P) D Ker(Jo). Together these two results prove that Ker(P) = Ker(Jo), so obviously 
the rank of the two matrices must be equal. Since rank(P) = N — m, this is also the rank 
of Jo, which completes the proof. 

In short, we have shown that in the limit a —>■ 0, the Jacobi-matrix J has N — m negative 
eigenvalues. 

4. Number of positive eigenvalues for a —00 

We now turn to the limit a —00. In this limit we rewrite J as follows: J = F'{n )aNL- 
D. Here D = diag((T~^, . . . , o"^^, 1, . . . ,1), which in the limit cr —00 becomes: 

lim D = diag(0, . . . , 0, 1, . . . , 1) = Q (29) 

Again, Q is a projection matrix, which now projects to the subspace spanned by the 
last m unit vectors, so Q is complementary to P. Following the same line of reasoning, but 
keeping in mind that now the constant factor in front of J_oo is negative, we find that in 
the limit a ^ — 00, the matrix J has m positive eigenvalues. 
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5. Number of positive eigenvalues of J for general a 

We are now ready to draw the conclusion. Just below o" = the matrix J must, by 
continuity, have at least N — m negative eigenvalues. If we now move from towards — oo, 
beyond a certain point there must be at least m positive eigenvalues (or equivalently, at 
most N — m — 1 negative eigenvalues). We already know [cf. Eq. (^)] that along the way 
exactly one eigenvalue changes sign, at a = —m/ {N — m). Taken together, this means that 
J has m positive and N — m — 1 negative eigenvalues for a < —m/{N — m), and m — 1 
positive and N — m negative eigenvalues for a > —m/{N — m). 

This completes the determination of the number of positive eigenvalues for the various 
branches in the bifurcation diagram. 
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